#### Replication file for "Ideological Moderation and Success in U.S. Elections, 2020-2022" 
#### Journal of Politics

#### Instructions
# 1) Load files 
#   a) Function files into subfolder called h_functions
#   b) Data files into subfolder called Data

# 2) Run code below

## Load packages
library(tidyverse)
library(estimatr)
library(tinytex)
library(gridExtra)
library(readxl)
library(scales)
library(colorspace)
library(kableExtra)
library(haven)
library(fixest)

## Load functions
sapply(paste0("h_functions/", list.files("h_functions")), source)


######################
## Figure 1  [No calculations]
######################

######################
## Table 1: Predicting legislative ideology 2022
######################

# Load Bonica scores
mgMCI = read.csv(file = "Data/mergeBonica.csv")

## Predict 2022 legislative ideology as a function of 
mci2020 <- read_excel("Data/main_candidate_information_2020Congress.xlsx", 
                      na= "NA") %>% 
  mutate(partyVAR = party, 
         "state" = substr(stateDist, 1, 2)) %>%
  filter(!is.na(nominateAdj)) %>%
  dplyr::select(c("candStem2020" = candStem, candidate, lastName, partyVAR, state, 
                  stateDist, nominateAdj))

mci2022 <- read_excel("Data/main_candidate_information_2022Congress.xlsx", na= "NA") %>% 
  mutate(partyVAR = party, 
         "state" = substr(stateDist, 1, 2)) %>%
  filter(!is.na(nominateAdj)) %>%
  dplyr::select(c("candStem2022" = candStem, candidate, lastName, partyVAR, state, 
                  stateDist, nominateAdj))

idealAll2020 = read_csv("Data/IdealPointsAll2020_labelFish.csv") %>%
  mutate("avgIdeal2020" = avgIdeal,
         "webIdeal2020" = idealWebLibcon, 
         "twIdeal2020"  = idealTwitterLibcon,
         "state" = substr(stateDist, 1, 2),
         "candStemLength" = nchar(candStem)) %>%
  mutate_at(c("avgIdeal2020", "twWeight", "thetaExternal"), round, 2) %>%
  mutate(candidateParty = paste0(candidate," (", party,")"),
         district = case_when(str_detect(district, "-Se") ~ paste0(district, "n"),
                              TRUE ~ district))   
# Merge in 2020 ideology 
mergeTEMP = mci2022 %>% 
  left_join(idealAll2020, by = "candidate") %>% 
  left_join(mgMCI %>% dplyr::select(c(candidate, recipient.cfscore)), by = "candidate") 
mergeTEMP$newMember = (!mci2022$candidate %in% mci2020$candidate)

# Republicans
all2022_avgR <- fixest::feols(scale(nominateAdj) ~ scale(avgIdeal2020),
                              data = mergeTEMP %>% filter(party == "R"))
all2022_bonR <- fixest::feols(scale(nominateAdj) ~ scale(recipient.cfscore),
                              data = mergeTEMP %>% filter(party == "R"))
new2022_avgR <- fixest::feols(scale(nominateAdj) ~ scale(avgIdeal2020),
                              data = mergeTEMP %>% filter(party == "R" & newMember ==1))
new2022_bonR <- fixest::feols(scale(nominateAdj) ~ scale(recipient.cfscore),
                              data = mergeTEMP %>% filter(party == "R" & newMember ==1))
# Democrats
all2022_avgD <- fixest::feols(scale(nominateAdj) ~ scale(avgIdeal2020),
                              data = mergeTEMP %>% filter(party == "D"))
all2022_bonD <- fixest::feols(scale(nominateAdj) ~ scale(recipient.cfscore),
                              data = mergeTEMP %>% filter(party == "D"))
new2022_avgD <- fixest::feols(scale(nominateAdj) ~ scale(avgIdeal2020),
                              data = mergeTEMP %>% filter(party == "D" & newMember ==1))
new2022_bonD <- fixest::feols(scale(nominateAdj) ~ scale(recipient.cfscore),
                              data = mergeTEMP %>% filter(party == "D" & newMember ==1))

# Set style to omit the "Fit statistics" label everywhere
sty <- style.tex(tablefoot.value = "Standard errors in parentheses; significance codes: *** p<0.01, ** p<0.05, * p<0.10.",
                 stats.title = "", var.title = "")

fixest::etable(list(all2022_bonR, all2022_avgR,
                    new2022_bonR, new2022_avgR,
                    all2022_bonD, all2022_avgD,
                    new2022_bonD, new2022_avgD), 
               tex=TRUE, digits = 2,
               fontsize = "footnotesize", 
               title = "Predicting legislative ideology, 2022", 
               label = "tab:repCompareDonor",  
               se.below = TRUE, #se.prefix = "[", se.suffix = "]",
               fitstat=c('n', 'r2'), digits.stats = 'r2',
               drop = c("Constant"),
               depvar = FALSE,
               order = "!Constant",
               style.tex = sty,
               dict = c("duckM" = "Legislative ideology 2020",
                        "scale(nominateAdj)" = "Legislative ideology 2022",
                        "scale(avgIdeal2020)" = "Text ideol. 2020",
                        "Bonicascores" = "Bonica scores",
                        "scale(recipient.cfscore)" = "CF ideol. 2020",
                        "Avg.textscores" = "Text scores 2020",
                        "Twittertextscores" = "Twitter text scores",
                        "Webtextscores" = "Web text scores",
                        "dem_inc" = "Dem. House inc."),
               notes = "", #\\hspace{5em}Standard errors in parentheses
               headers = list(
                 # Top row: party
                 ":_:" = list("Republicans" = 4, "Democrats" = 4),
                 #":_:Party" = list("Republicans" = 4, "Democrats" = 4),
                 # Second row: within-party groups
                 " " = list("All" = 2, "New members" = 2,
                            "All" = 2, "New members" = 2)  ) )

######################
## Figure 2  [No calculations]
######################

######################
## Figure 3: Ideology and Democratic House performance compared to Biden
######################
r_hp_2020 = kappa_eachYear_House_FUNCTION(2020, compare = "Pres", cluster=1)
r_hp_2022 = kappa_eachYear_House_FUNCTION(2022, compare = "Pres", cluster=1)
r_hs_2020 = kappa_eachYear_House_FUNCTION(2020, compare = "Sen",  cluster=1)

yearList = c(2020, 2022, 2020, 2022)

par(mfrow=c(2, 2), 
    mar=c(1.0, 2.0, 0.75,  0.2), 
    oma=c(1.1, 1.8, 0.5,  0.2))	 	# mar(south, west, north, east)
for(plotCount in 1:4){
  if(plotCount==4){
    plot(c(-0.6, 0.3), c(-10, 8), 
         xlab = "", ylab = "", 
         type="n", cex.axis=0.9, 
         xaxt="n", yaxt="n", ylim = c(-11.4, 8.4))
    axis(1, tick = T, tck = -0.03, 
         at = seq(-1, 1, by=0.25), 
         labels = seq(-1, 1, by=0.25),  
         cex.axis = 0.7, mgp = c(2, 0.0, 0))
    axis(2, las = 1, tick = T, tck = -0.03, 
         cex.axis = 0.7, 
         mgp = c(2, 0.25, 0), 
         at = seq(-10, 10, by=2), 
         labels = paste0(seq(-10, 10, by=2), "%")) 
    text(-0.14, -1, "Data unavailable", col = "darkgrey")}
  
  if(plotCount!=4){
    if(plotCount == 1) r.tmp = r_hp_2020
    if(plotCount == 2) r.tmp = r_hp_2022
    if(plotCount == 3) r.tmp = r_hs_2020
    plot(r.tmp$mdiHouseAll$demDiff_HousePres ~r.tmp$mdiHouseAll$kappaDiff_HousePres, 
         xlab = "", ylab = "",
         type="n", cex.axis=0.9, 
         xaxt="n", yaxt="n", ylim = c(-11.4, 8.4))  #xlim = c(-0.5, 0.75)
    
    abline(r.tmp$all_no_covar,    col="darkgrey", lwd = 2)
    abline(r.tmp$compet_no_covar, col="#3c78c3",  lwd = 3)
    points(r.tmp$mdiHouseAll$demDiff_HousePres ~ r.tmp$mdiHouseAll$kappaDiff_HousePres, 
           pch=16, cex = 0.8, col="darkgrey")
    points(r.tmp$mdiHouseCompet$demDiff_HousePres ~ r.tmp$mdiHouseCompet$kappaDiff_HousePres,
           pch=16, cex = 1.0, col="#3c78c3")
    
    if(plotCount< 3) {
      title(yearList[plotCount], line = 0.19, cex.main =0.9)}
    
    if(plotCount ==3){
      mtext("House\nDem.\nmargin\nover\nSenate\nDem.", side = 2, 
            las = 1, at = 0, line = 1.2, cex = 0.75)
      mtext("Cutpoint difference", side = 1, 
            las = 1, 
            at = 0.72, 
            line = 1.1, cex = 1)}
    # KEY
    if(plotCount ==1){
      points(0.35, -10.0, pch=16, cex = 1.0, col="#3c78c3")
      text  (0.48, -10.0, "Competitive", col="#3c78c3", cex = 0.7)
      points(0.35, -11.2, pch=16, cex = 0.8, col="darkgrey")
      text  (0.51, -11.2, "Not competitive", col="darkgrey", cex = 0.7) }
    if(plotCount ==2){
      points(-0.60, 8.3, pch=16, cex = 1.0, col="#3c78c3")
      text  (-0.49, 8.3, "Competitive", col="#3c78c3", cex = 0.7)
      points(-0.60, 7, pch=16, cex = 0.8, col="darkgrey")
      text  (-0.46, 7, "Not competitive", col="darkgrey", cex = 0.7) }
    if(plotCount ==3){
      points(0.35, -10.0, pch=16, cex = 1.0, col="#3c78c3")
      text  (0.47, -10.0, "Competitive", col="#3c78c3", cex = 0.7)
      points(0.35, -11.2, pch=16, cex = 0.8, col="darkgrey")
      text  (0.50, -11.2, "Not competitive", col="darkgrey", cex = 0.7) }
    
    axis(1, tick = T, tck = -0.03, 
         at = seq(-1, 1, by=0.25), 
         labels = seq(-1, 1, by=0.25),  
         cex.axis = 0.5, mgp = c(2, 0.0, 0))
    axis(2, las = 1, tick = T, tck = -0.03, 
         cex.axis = 0.5, 
         mgp = c(2, 0.25, 0), 
         at = seq(-10, 10, by=2), 
         labels = paste0(seq(-10, 10, by=2), "%")) 
    if(plotCount == 1){
      mtext("House\nDem.\nmargin\nover\nBiden\n2020", side = 2, 
            las = 1, at = 0, line = 1.2, cex = 0.75)}
  } }


######################
## Figure 4: Ideology and Elections: Comparisons to President and Senate
######################

par(mfrow=c(1, 2),
    mar=c(1.0, 1.2, 3, 0.5),
    oma=c(1.0, 4,   2, 0.3))  # mar(south, west, north, east)
betaFigure_House_FUNCTION(year = 2020)
betaFigure_House_FUNCTION(year = 2022)


######################
## Figure 5: Counterfactuals
######################
par(mfrow=c(1, 1),
    mar=c(1.0, 1.2, 3, 0.5),
    oma=c(1.0, 4,   2, 0.4))  # mar(south, west, north, east)

## Create figure for counterfactual
cf2020 = simChange_FUNCTION(2020) %>% unlist()
cf2022 = simChange_FUNCTION(2022) %>% unlist()

cf = rbind(c(cf2020[1:3], cf2020[1], cf2020[4:5]), 
           c(cf2022[1:3], cf2022[1], cf2022[4:5]))

plot(1:6, seq(min(cf)-8, max(cf)+20, length = 6),
     type="n", lty=2, xlab="", ylab="", xaxt='n', cex.axis=0.9, yaxt='n')
cf.labels = c("Actual\n", "Dem. moderate\nRep. actual", "Dem. moderate\nRep. extreme",
              "Actual\n", "Dem. actual\nRep. moderate", "Dem. extreme\nRep. moderate")
axis(1, at=seq(1:6),	labels=cf.labels, cex.axis=0.6, padj = -0.5,  tck = -0.02)
axis(2, at=seq(180, 250, by = 10),	labels=seq(180, 250, by = 10), 
     cex.axis=0.65, las = T,   mgp = c(3, 0.5, 0), tck = -0.025)
abline(v= 3.5, col="grey")
abline(h= 218, col="grey")
text(2.8, 214, "218 seats required for\nHouse majority", col= "grey", cex=0.6)
mtext("Number",    side = 2, las= T, cex = 0.8, line = 2.2, at = 240)
mtext("of",        side = 2, las= T, cex = 0.8, line = 3.0, at = 234)
mtext("Democrats", side = 2, las= T, cex = 0.8, line = 1.8, at = 228)


# Set colors
blues <- colorRampPalette(c("lightblue", "blue", "darkblue"))
#barplot(1:10, col = blues(10), main = "Blue Color Scale")
color_dem = blues(10)[c(4, 10)]
reds <- colorRampPalette(c("red", "tomato2", "darkred"))
color_rep = reds(10) [c(4, 10)]

for(rr in 1:2){
  lines(1:3, cf[rr,1:3], col= color_dem[rr])
  points(1:3, cf[rr,1:3], pch = 19, col= color_dem[rr])
  
  lines(4:6, cf[rr,4:6], col= color_rep[rr])
  points(4:6, cf[rr,4:6], pch = 19, col= color_rep[rr])
}

text(2.1, 255, "Changes Helping Democrats", col= "darkblue", cex = 0.8)
text(4.9, 255, "Changes Helping Republicans", col= "darkred", cex = 0.8)

text(1.2 , 225, "2020", col= color_dem[1], cex = 0.7)
text(1.23, 212, "2022", col= color_dem[2], cex = 0.7)

text(4.2 , 224, "2020", col= color_rep[1], cex = 0.7)
text(4.13, 210, "2022", col= color_rep[2], cex = 0.7)

######################
## Figure 6: Ideology and Democratic Senate, presidential and gubernatorial races
######################
# Senate-President comparisons
county_res_sp = kappa_bothYears_SenPresGov_FUNCTION(geography = "County", compare = "sp")
state_res_sp  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "sp")

# Senate-Governor comparisons
county_res_sg = kappa_bothYears_SenPresGov_FUNCTION(geography = "County", compare = "sg")
state_res_sg  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "sg")

# President - Governor comparisons
county_res_pg = kappa_bothYears_SenPresGov_FUNCTION(geography = "County", compare = "pg")
state_res_pg  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "pg")

percentile <- function(x) round((rank(x) - 1)/(length(x) - 1), 2)
par(mfrow=c(3, 2),
    mar=c(1, 1.2, 1.,0.5),
    oma=c(2.0, 4.0,  1, 0.3))  # mar(south, west, north, east)
for(plotCount in 1:6){
  if(plotCount ==1) {
    tmp = county_res_sp
    yMin = -12
    yMax =  12}
  if(plotCount ==2) {
    tmp = state_res_sp
    yMin = -12
    yMax =  12}
  if(plotCount ==3) {
    tmp = county_res_sg
    yMin = -25
    yMax =  25}
  if(plotCount ==4) {
    tmp = state_res_sg
    yMin = -25
    yMax =  25}
  if(plotCount ==5) {
    tmp = county_res_pg
    yMin = -20
    yMax =  50}
  if(plotCount ==6) {
    tmp = state_res_pg
    yMin = -20
    yMax =  50}
  
  plot(100*tmp$dta$dem_diff ~ tmp$dta$cutpoint_diff, 
       xlab = "", ylab = "",
       type="n", cex.axis=0.9, 
       xaxt="n", yaxt="n", ylim = c(yMin, yMax))  #, ylim = c(-56, 12), xlim = c(-0.5, 0.75)
  #c(min(tmp$dta$dem_diff, na.rm= T), max(tmp$dta$dem_diff, na.rm= T))
  if(plotCount == 1) mtext("Senate\nDem.\nmargin\nover\nBiden", side = 2,
                           las = 1, at = 0, line = 1.6, cex = 0.75)
  if(plotCount == 3) mtext("Senate\nDem.\nmargin\nover\nDem.\nGov.", side = 2, 
                           las = 1, at = 0, line = 1.6, cex = 0.75)
  if(plotCount == 5) mtext("Biden\n2020\nmargin\nover\nDem.\nGov.", side = 2, 
                           las = 1, at = 15, line = 1.6, cex = 0.75)
  if(plotCount==1) title("County", line = 0.3, cex.main =1)
  if(plotCount==2) title("State",  line = 0.3, cex.main =1)
  
  abline(tmp$res_no_covar, col="#3c78c3", lwd = 2)
  if(plotCount %in% c(1,3,5)) points(100*tmp$dta$dem_diff ~ tmp$dta$cutpoint_diff, 
                                     pch=16, cex = percentile(tmp$dta$county_votes)*1.4, col="#3c78c3")
  if(plotCount %in% c(2,4,6)) points(100*tmp$dta$dem_diff ~ tmp$dta$cutpoint_diff, 
                                     pch=16, cex = 1.2, col="#3c78c3")
  axis(1, tick = T, tck = -0.03, 
       at = seq(-1, 1, by=0.25), 
       labels = seq(-1, 1, by=0.25),  
       cex.axis = 0.5, mgp = c(2, 0.0, 0))
  axis(2, las = 1, tick = T, tck = -0.03, 
       cex.axis = 0.5, 
       mgp = c(2, 0.25, 0), 
       at = seq(-70, 70, by=10), 
       labels = paste0(seq(-70, 70, by=10), "%")) 
  if(plotCount> 4) mtext("Cutpoint difference", side = 1, 
                         las = 1, at = 0.2, line = 1.4, cex = 0.8)
}


######################
## Figure 7: Ideology and elections: Comparisons of Senate, presidential and gub. results
######################

par(mfrow=c(1, 1),
    mar=c(1.0, 2.2, 3, 1.5),
    oma=c(2.0, 4,   2, 1.3))  # mar(south, west, north, east)

betaFigure_SenPresGov_FUNCTION()

######################
## Table 2: Ideology and incumbency interactions
######################
# Ideology x incumbency: House Pres in 2020
allTEMP = r_hp_2020$mdiHouseAll
compTEMP= r_hp_2020$mdiHouseCompet

formula_incDiff  <- as.formula("demDiff_HousePres ~ dem_inc + rep_inc+ 
                               dem_avgIdeal*dem_inc + rep_avgIdeal*rep_inc")

incDiff20       <- fixest::feols(formula_incDiff,
                                 data = allTEMP, cluster ~ state)
#se = "iid")
summary(incDiff20)
incDiff_c20       <- fixest::feols(formula_incDiff,
                                   data = compTEMP, cluster ~ state) # ,se = "iid")
summary(incDiff_c20)

# Ideology x incumbency: House Pres in 2022
allTEMP= r_hp_2022$mdiHouseAll
compTEMP= r_hp_2022$mdiHouseCompet
formula_incDiff  <- as.formula("demDiff_HousePres ~ dem_inc + rep_inc+ dem_avgIdeal*dem_inc + rep_avgIdeal*rep_inc")
incDiff22       <- fixest::feols(formula_incDiff,       
                                 data = allTEMP, cluster ~ state) # se = "iid")
summary(incDiff22)

incDiff_c22       <- fixest::feols(formula_incDiff,       
                                   data = compTEMP, cluster ~ state) # se = "iid")
summary(incDiff_c22)

# Ideology x incumbency: Senate Pres in 2020-2022
spTEMP  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "sp", weights=1)

formula_incDiff_SenPres  <- as.formula("dem_diff*100 ~ dem_sen_avgIdeal*dem_sen_inc +
rep_sen_avgIdeal*rep_sen_inc + dem_sen_inc + rep_sen_inc + `2020`")
thetaDiff_all       <- fixest::feols(formula_incDiff_SenPres,       
                                     data = spTEMP$dta, se = "iid") 
summary(thetaDiff_all)

fixest::etable(list(incDiff20, incDiff_c20, incDiff22, incDiff_c22, thetaDiff_all),
               tex=TRUE, digits = 2,
               fontsize = "footnotesize", 
               title = "Ideology and incumbency interactions", 
               label = "tab:incumbencyInteractions",  
               se.below = F, fitstat=c('n', 'r2'), digits.stats = 'r2',
               order = "!Constant",
               drop = c("Constant", '2020'),
               style.tex = sty,
               dict = c("dem_inc" = "Dem. incumbent",
                        "rep_inc" = "Rep. incumbent",
                        "dem_avgIdeal" = "Dem. ideology", 
                        "rep_avgIdeal" = "Rep. ideology", 
                        "dem_inc:dem_avgIdeal" = "Dem. ideol. x incumbent", 
                        "rep_inc:rep_avgIdeal" = "Rep. ideol. x incumbent", 
                        "dem_sen_avgIdeal"     = "Dem. ideology", 
                        "dem_sen_inc"          = "Dem. incumbent",
                        "rep_sen_avgIdeal"     = "Rep. ideology", 
                        "rep_sen_inc"          = "Rep. incumbent",
                        "'2020'"               = "2020", 
                        "dem_sen_avgIdeal:dem_sen_inc" = "Dem. ideol. x incumbent", 
                        "rep_sen_avgIdeal:rep_sen_inc" = "Rep. ideol. x incumbent", 
                        "dem_diff * 100"    = "Difference in Democratic percent",
                        "demDiff_HousePres" = "Difference in Democratic percent"),
               notes = "",
               headers = list(":_:" = list("House-Pres. 2020" = 2, 
                                           "House-Pres. 2022" = 2, "Sen.-Pres." = 1)))

######################
## Table 3: Party Differences Cutpoint Comparisons
######################

# Assess whether coefficient varies by party: House-President 2020
r_hp_2020 = kappa_eachYear_House_FUNCTION(2020, compare = "Pres", measure = "avgText", clusterCode =1)

all = r_hp_2020$mdiHouseAll %>%
  mutate(thetaDiff_DemHousePres = (dem_avgIdeal - dem_pres) / 2,
         thetaDiff_RepHousePres = (rep_avgIdeal - rep_pres) / 2)
compet = r_hp_2020$mdiHouseCompet %>%
  mutate(thetaDiff_DemHousePres = (dem_avgIdeal - dem_pres) / 2,
         thetaDiff_RepHousePres = (rep_avgIdeal - rep_pres) / 2)

formula_thetaDiff_inc_covar  <- as.formula("demDiff_HousePres ~ thetaDiff_DemHousePres + thetaDiff_RepHousePres + 
                                           house_dem_inc + house_rep_inc")
thetaDiff_hp20_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, cluster ~ state )
summary(thetaDiff_hp20_all)

thetaDiff_hp20_compet       <- fixest::feols(formula_thetaDiff_inc_covar,       data = compet, cluster ~ state)
summary(thetaDiff_hp20_compet)

# Assess whether coefficient varies by party: House-President 2022
r_hp_2022 = kappa_eachYear_House_FUNCTION(2022, compare = "Pres", clusterCode =1)
names(r_hp_2022$mdiHouseAll)[grep("diff", names(r_hp_2020$mdiHouseAll), ignore.case = TRUE)]

all = r_hp_2022$mdiHouseAll %>%
  mutate(thetaDiff_DemHousePres = (dem_avgIdeal - dem_pres) / 2,
         thetaDiff_RepHousePres = (rep_avgIdeal - rep_pres) / 2)
compet = r_hp_2022$mdiHouseCompet %>%
  mutate(thetaDiff_DemHousePres = (dem_avgIdeal - dem_pres) / 2,
         thetaDiff_RepHousePres = (rep_avgIdeal - rep_pres) / 2)

formula_thetaDiff_inc_covar  <- as.formula("demDiff_HousePres ~ thetaDiff_DemHousePres + thetaDiff_RepHousePres + 
                                           house_dem_inc + house_rep_inc")
thetaDiff_hp22_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, cluster ~ state )
summary(thetaDiff_hp22_all)

thetaDiff_hp22_compet       <- fixest::feols(formula_thetaDiff_inc_covar,       data = compet, cluster ~ state)
summary(thetaDiff_hp22_compet)

# Assess whether coefficient varies by party: House-Senate 2020
r_hs_2020 = kappa_eachYear_House_FUNCTION(2020, compare = "Sen",  clusterCode =1)

all = r_hs_2020$mdiHouseAll %>%
  mutate(thetaDiff_DemHouseSen = (dem_house_avgIdeal - dem_sen_avgIdeal) / 2,
         thetaDiff_RepHouseSen = (rep_house_avgIdeal - rep_sen_avgIdeal) / 2, 
         dem_sen_inc = 1*(dem_sen_inc==1),
         rep_sen_inc = 1*(rep_sen_inc==1))

formula_thetaDiff_inc_covar  <- as.formula("demDiff_HouseSen ~ thetaDiff_DemHouseSen + thetaDiff_RepHouseSen + 
                                           house_dem_inc + house_rep_inc + dem_sen_inc + rep_sen_inc")
thetaDiff_hs20_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, cluster ~ state )
summary(thetaDiff_hs20_all)

compet = r_hs_2020$mdiHouseCompet %>%
  mutate(thetaDiff_DemHouseSen = (dem_house_avgIdeal - dem_sen_avgIdeal) / 2,
         thetaDiff_RepHouseSen = (rep_house_avgIdeal - rep_sen_avgIdeal) / 2, 
         dem_sen_inc = 1*(dem_sen_inc==1),
         rep_sen_inc = 1*(rep_sen_inc==1))


thetaDiff_hs20_compet       <- fixest::feols(formula_thetaDiff_inc_covar,       data = compet, cluster ~ state)
summary(thetaDiff_hs20_compet)

# Assess whether coefficient varies by party: Senate Pres 2020-2022
state_res_sp  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "sp", weights=1)

names(state_res_sp$dta)[grep("020", names(state_res_sp$dta), ignore.case = TRUE)]

all = state_res_sp$dta %>%
  mutate(thetaDiff_DemSenPres = (dem_sen_avgIdeal) / 2,
         thetaDiff_RepSenPres = (rep_sen_avgIdeal) / 2, 
         dem_sen_inc = 1*(dem_sen_inc==1),
         rep_sen_inc = 1*(rep_sen_inc==1),
         dem_diff = 100*dem_diff)

formula_thetaDiff_inc_covar  <- as.formula("dem_diff ~ thetaDiff_DemSenPres + thetaDiff_RepSenPres + 
                                           dem_sen_inc + rep_sen_inc + `2020`")
thetaDiff_sp_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, cluster ~ state )
summary(thetaDiff_sp_all)

# Assess whether coefficient varies by party: Senate Gov 2020-2022
state_res_sg  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "sg", weights=1)

names(state_res_sg$dta)[grep("diff", names(state_res_sg$dta), ignore.case = TRUE)]

all = state_res_sg$dta %>%
  mutate(thetaDiff_DemSenGov = (dem_sen_avgIdeal - dem_gov_avgIdeal) / 2,
         thetaDiff_RepSenGov = (rep_sen_avgIdeal - rep_gov_avgIdeal) / 2, 
         dem_sen_inc = 1*(dem_sen_inc==1),
         rep_sen_inc = 1*(rep_sen_inc==1),
         dem_diff = 100*dem_diff)

formula_thetaDiff_inc_covar  <- as.formula("dem_diff ~ thetaDiff_DemSenGov + thetaDiff_RepSenGov + 
                                           dem_sen_inc + rep_sen_inc + dem_gov_inc + rep_gov_inc + `2020`")
thetaDiff_sg_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, weights = ~ county_votes)
summary(thetaDiff_sg_all)

# Assess whether coefficient varies by party: Senate Gov 2020-2022
state_res_pg  = kappa_bothYears_SenPresGov_FUNCTION(geography = "State",  compare = "pg", weights=1)

names(state_res_pg$dta)[grep("pres", names(state_res_pg$dta), ignore.case = TRUE)]

all = state_res_pg$dta %>%
  mutate(thetaDiff_DemPresGov = (dem_pres-dem_gov_avgIdeal) / 2,
         thetaDiff_RepPresGov = (rep_pres-rep_gov_avgIdeal) / 2, 
         dem_diff = 100*dem_diff)
formula_thetaDiff_inc_covar  <- as.formula("dem_diff ~ thetaDiff_DemPresGov + thetaDiff_RepPresGov + 
                                           dem_gov_inc + rep_gov_inc + `2020`")
thetaDiff_pg_all       <- fixest::feols(formula_thetaDiff_inc_covar,       data = all, weights = ~ county_votes)
summary(thetaDiff_pg_all)

## Create a data frame that contains information for the table
partyDiff_table = data.frame(DemCut= rep(NA, 6), DemCut_t= rep(NA, 6), 
                             RepCut= rep(NA, 6), RepCut_t= rep(NA, 6),
                             DemCutCompet= rep(NA, 6), DemCutCompet_t= rep(NA, 6), 
                             RepCutCompet= rep(NA, 6), RepCutCompet_t= rep(NA, 6))
row.names(partyDiff_table) = c("HousePres2020", "HousePres2022", "HouseSen2020",
                           "SenPres", "SenGov", "PresGov")

partyDiff_table["HousePres2020", "DemCut"] =
  round(summary(thetaDiff_hp20_all)$coeftable["thetaDiff_DemHousePres", "Estimate"], 2)
partyDiff_table["HousePres2020", "DemCut_t"] =
  round(summary(thetaDiff_hp20_all)$coeftable["thetaDiff_DemHousePres", "t value"], 2)
partyDiff_table["HousePres2020", "RepCut"] =
  round(summary(thetaDiff_hp20_all)$coeftable["thetaDiff_RepHousePres", "Estimate"], 2)
partyDiff_table["HousePres2020", "RepCut_t"] =
  round(summary(thetaDiff_hp20_all)$coeftable["thetaDiff_RepHousePres", "t value"], 2)

partyDiff_table["HousePres2020", "DemCutCompet"] =
  round(summary(thetaDiff_hp20_compet)$coeftable["thetaDiff_DemHousePres", "Estimate"], 2)
partyDiff_table["HousePres2020", "DemCutCompet_t"] =
  round(summary(thetaDiff_hp20_compet)$coeftable["thetaDiff_DemHousePres", "t value"], 2)
partyDiff_table["HousePres2020", "RepCutCompet"] =
  round(summary(thetaDiff_hp20_compet)$coeftable["thetaDiff_RepHousePres", "Estimate"], 2)
partyDiff_table["HousePres2020", "RepCutCompet_t"] =
  round(summary(thetaDiff_hp20_compet)$coeftable["thetaDiff_RepHousePres", "t value"], 2)

partyDiff_table["HousePres2022", "DemCut"] =
  round(summary(thetaDiff_hp22_all)$coeftable["thetaDiff_DemHousePres", "Estimate"], 2)
partyDiff_table["HousePres2022", "DemCut_t"] =
  round(summary(thetaDiff_hp22_all)$coeftable["thetaDiff_DemHousePres", "t value"], 2)
partyDiff_table["HousePres2022", "RepCut"] =
  round(summary(thetaDiff_hp22_all)$coeftable["thetaDiff_RepHousePres", "Estimate"], 2)
partyDiff_table["HousePres2022", "RepCut_t"] =
  round(summary(thetaDiff_hp22_all)$coeftable["thetaDiff_RepHousePres", "t value"], 2)

partyDiff_table["HousePres2022", "DemCutCompet"] =
  round(summary(thetaDiff_hp22_compet)$coeftable["thetaDiff_DemHousePres", "Estimate"], 2)
partyDiff_table["HousePres2022", "DemCutCompet_t"] =
  round(summary(thetaDiff_hp22_compet)$coeftable["thetaDiff_DemHousePres", "t value"], 2)
partyDiff_table["HousePres2022", "RepCutCompet"] =
  round(summary(thetaDiff_hp22_compet)$coeftable["thetaDiff_RepHousePres", "Estimate"], 2)
partyDiff_table["HousePres2022", "RepCutCompet_t"] =
  round(summary(thetaDiff_hp22_compet)$coeftable["thetaDiff_RepHousePres", "t value"], 2)

partyDiff_table["HouseSen2020", "DemCut"] =
  round(summary(thetaDiff_hs20_all)$coeftable["thetaDiff_DemHouseSen", "Estimate"], 2)
partyDiff_table["HouseSen2020", "DemCut_t"] =
  round(summary(thetaDiff_hs20_all)$coeftable["thetaDiff_DemHouseSen", "t value"], 2)
partyDiff_table["HouseSen2020", "RepCut"] =
  round(summary(thetaDiff_hs20_all)$coeftable["thetaDiff_RepHouseSen", "Estimate"], 2)
partyDiff_table["HouseSen2020", "RepCut_t"] =
  round(summary(thetaDiff_hs20_all)$coeftable["thetaDiff_RepHouseSen", "t value"], 2)

partyDiff_table["HouseSen2020", "DemCutCompet"] =
  round(summary(thetaDiff_hs20_compet)$coeftable["thetaDiff_DemHouseSen", "Estimate"], 2)
partyDiff_table["HouseSen2020", "DemCutCompet_t"] =
  round(summary(thetaDiff_hs20_compet)$coeftable["thetaDiff_DemHouseSen", "t value"], 2)
partyDiff_table["HouseSen2020", "RepCutCompet"] =
  round(summary(thetaDiff_hs20_compet)$coeftable["thetaDiff_RepHouseSen", "Estimate"], 2)
partyDiff_table["HouseSen2020", "RepCutCompet_t"] =
  round(summary(thetaDiff_hs20_compet)$coeftable["thetaDiff_RepHouseSen", "t value"], 2)

partyDiff_table["SenPres", "DemCut"] =
  round(summary(thetaDiff_sp_all)$coeftable["thetaDiff_DemSenPres", "Estimate"], 2)
partyDiff_table["SenPres", "DemCut_t"] =
  round(summary(thetaDiff_sp_all)$coeftable["thetaDiff_DemSenPres", "t value"], 2)
partyDiff_table["SenPres", "RepCut"] =
  round(summary(thetaDiff_sp_all)$coeftable["thetaDiff_RepSenPres", "Estimate"], 2)
partyDiff_table["SenPres", "RepCut_t"] =
  round(summary(thetaDiff_sp_all)$coeftable["thetaDiff_RepSenPres", "t value"], 2)

partyDiff_table["SenGov", "DemCut"] =
  round(summary(thetaDiff_sg_all)$coeftable["thetaDiff_DemSenGov", "Estimate"], 2)
partyDiff_table["SenGov", "DemCut_t"] =
  round(summary(thetaDiff_sg_all)$coeftable["thetaDiff_DemSenGov", "t value"], 2)
partyDiff_table["SenGov", "RepCut"] =
  round(summary(thetaDiff_sg_all)$coeftable["thetaDiff_RepSenGov", "Estimate"], 2)
partyDiff_table["SenGov", "RepCut_t"] =
  round(summary(thetaDiff_sg_all)$coeftable["thetaDiff_RepSenGov", "t value"], 2)

partyDiff_table["PresGov", "DemCut"] =
  round(summary(thetaDiff_pg_all)$coeftable["thetaDiff_DemPresGov", "Estimate"], 2)
partyDiff_table["PresGov", "DemCut_t"] =
  round(summary(thetaDiff_pg_all)$coeftable["thetaDiff_DemPresGov", "t value"], 2)
partyDiff_table["PresGov", "RepCut"] =
  round(summary(thetaDiff_pg_all)$coeftable["thetaDiff_RepPresGov", "Estimate"], 2)
partyDiff_table["PresGov", "RepCut_t"] =
  round(summary(thetaDiff_pg_all)$coeftable["thetaDiff_RepPresGov", "t value"], 2)

partyDiff_table

######################
## Figure 8: Effect of Ideology Over Time
######################

YearList= seq(2004, 2020, by = 4)
res_dime_year= data.frame("Year" = YearList, 
                          "dwdime" = NA, "dwdime_uci" = NA, "dwdime_lci" = NA, 
                          "dwdime_n" = NA,
                          "cf" = NA, "cf_uci" = NA, "cf_lci" = NA, "cf_n" = NA)
for(yy in YearList){
  dw_temp = readRDS(paste0("Data//dwDimeOLS_", yy, ".rds"))
  cf_temp = readRDS(paste0("Data//cfDimeOLS_", yy, ".rds"))
  
  res_dime_year[res_dime_year$Year == yy, "dwdime"] = dw_temp$coefficients["kappa_house"]
  res_dime_year[res_dime_year$Year == yy, "dwdime_uci"] = dw_temp$coefficients["kappa_house"] + 
    1.96*summary(dw_temp)$coefficients["kappa_house", "Std. Error"]
  res_dime_year[res_dime_year$Year == yy, "dwdime_lci"] = dw_temp$coefficients["kappa_house"] - 
    1.96*summary(dw_temp)$coefficients["kappa_house", "Std. Error"]
  res_dime_year[res_dime_year$Year == yy, "dwdime_n"]= sum(summary(dw_temp)$df[1:2])
  
  res_dime_year[res_dime_year$Year == yy, "cf"]     = cf_temp$coefficients["kappa_house_cf"]
  res_dime_year[res_dime_year$Year == yy, "cf_uci"] = cf_temp$coefficients["kappa_house_cf"] + 
    1.96*summary(cf_temp)$coefficients["kappa_house_cf", "Std. Error"]
  res_dime_year[res_dime_year$Year == yy, "cf_lci"] = cf_temp$coefficients["kappa_house_cf"] - 
    1.96*summary(cf_temp)$coefficients["kappa_house_cf", "Std. Error"]
  res_dime_year[res_dime_year$Year == yy, "cf_n"]= sum(summary(cf_temp)$df[1:2])
}
res_dime_year

par(mfrow=c(1, 1),
    mar=c(1.0, 2.0, 1.,0.5),
    oma=c(1.5, 4.0,  1, 0.3))  # mar(south, west, north, east)
plot(YearList, res_dime_year$dwdime, type="n", lty=2, xlab="", ylab="", xaxt='n',
     cex.axis=0.9, yaxt='n', ylim = c(0, 18.1))# xlim= c(-3.7, 3.7)
mtext("Year", side = 1, at = mean(YearList), line = 1, cex = 1)
axis(1, at=YearList,	labels=YearList, cex.axis=0.7, padj = -1.9, tick=TRUE, tck = -0.03)
axis(2, at=seq(0, 20, by= 2.5),	labels=seq(0, 20, by= 2.5), mgp = c(3, 0.55, 0),
     cex.axis=0.6, las = T)

mtext("Coefficient", side = 2, at = 13.75,   line = 1,  cex = 1.0, las = T)


ok   <- with(res_dime_year, complete.cases(dwdime, dwdime_lci, dwdime_uci)) & !is.na(YearList)
x    <- YearList[ok]
mu   <- res_dime_year$dwdime[ok]
lci  <- res_dime_year$dwdime_lci[ok]
uci  <- res_dime_year$dwdime_uci[ok]
ord  <- order(x)
x    <- x[ord]; mu <- mu[ord]; lci <- lci[ord]; uci <- uci[ord]

shade_col <- adjustcolor("lightblue", alpha.f = 0.15)  # transparency
polygon(
  x = c(x, rev(x)),
  y = c(uci, rev(lci)),
  col = shade_col, border = NA
)

ok_cf   <- with(res_dime_year, complete.cases(cf, cf_lci, cf_uci)) & !is.na(YearList)
x_cf    <- YearList[ok]
mu_cf   <- res_dime_year$cf[ok]
lci_cf  <- res_dime_year$cf_lci[ok]
uci_cf  <- res_dime_year$cf_uci[ok]
ord_cf  <- order(x_cf)
x_cf    <- x[ord_cf]; mu_cf <- mu[ord_cf]; lci_cf <- lci_cf[ord_cf]; uci_cf <- uci_cf[ord_cf]

shade_col <- adjustcolor("lightgreen", alpha.f = 0.15)  # transparency
polygon(
  x = c(x_cf, rev(x_cf)),
  y = c(uci_cf, rev(lci_cf)),
  col = shade_col, border = NA
)

lines(YearList, res_dime_year$dwdime, 
      lwd = 2, col = "darkblue", pch = 19)
lines(YearList, res_dime_year$dwdime_uci,
      lwd = 1, lty=1, col = "lightblue", pch = 19)
lines(YearList, res_dime_year$dwdime_lci,
      lwd = 1, lty=1, col = "lightblue", pch = 19)
points(YearList, res_dime_year$dwdime, 
       lwd = 2, col = "darkblue", pch = 19)
text(2009.6, 13.1, "DW-Dime", col= "darkblue")

lines(YearList, res_dime_year$cf, 
      lwd = 2, col = "darkgreen", pch = 19)
points(YearList, res_dime_year$cf, 
       lwd = 2, col = "darkgreen", pch = 23)
lines(YearList, res_dime_year$cf_uci, 
      lwd = 1, lty=2,col = "lightgreen", pch = 19)
lines(YearList, res_dime_year$cf_lci, 
      lwd = 1, lty=2, col = "lightgreen", pch = 19)
text(2005.8, 10.4, "CF", col= "darkgreen")
